#include <omp.h>
#include "stdlib.h"
#include "stdio.h"
#include "fftw3.h"
#include "calc_type.h"
#include "math_statistics_func.h"


void fftw3_transform(uint32_t n, float *in, fftwf_complex *out)
{
  fftwf_plan plan_forward;
  
  plan_forward = fftwf_plan_dft_r2c_1d(n, in, out, FFTW_ESTIMATE); 	//out正负频率
  fftwf_execute(plan_forward);
  
  fftwf_destroy_plan(plan_forward);
}

void fftw3_inv_transform(uint32_t n, fftwf_complex *in, float *out)
{
  fftwf_plan plan_backward;

  plan_backward = fftwf_plan_dft_c2r_1d(n, in, out, FFTW_ESTIMATE);
  fftwf_execute(plan_backward);

  fftwf_destroy_plan(plan_backward);
}

void fftw3_transform_out(uint32_t len, float *in, float *out)
{
	uint32_t nc;

    fftwf_complex fft_out[DATA_CALC_MAX_SIZE + 1] = {{0.0, 0.0}, {0.0, 0.0}};
    float amplify_out[DATA_CALC_MAX_SIZE] = {0.0, 0.0};

	nc = (len / 2) + 1;
	
	//#pragma omp critical
	fftw3_transform(len, in, fft_out);
	arm_cmplx_mag_f32((float32_t *)fft_out, amplify_out, nc - 1);

	out[0] = 0.0;
	for(uint32_t i = 1; i < (len / 2); i++)
      	out[i] = (2 * amplify_out[i]) / (float)(len);
}



void fftw3_transform_test(uint32_t len, float *in, float *out)
{
    fftwf_complex fft_out[(DATA_CALC_MAX_SIZE / 2) + 1] = {{0.0, 0.0}, {0.0, 0.0}};
    float amplify_out[DATA_CALC_MAX_SIZE / 2] = {0.0, 0.0};
	
	fftw3_transform(len, in, fft_out);
	fftw3_inv_transform(len, fft_out, amplify_out);

	for(uint32_t i = 0; i < len; i++)	//正变换，反变换，得到原始波形，这里除以的是N
      	out[i] = amplify_out[i] / (float)(len);
}


